Goals

Data Structures

• Review data structures section in workbook
• Cover vectors, matrices, arrays, lists, data frames, etc.
• Cover data frame manipulation using tidyverse functions and the pipe operator (%>% or |>)

Why R for Spatial Data Analysis?

• Lightweight, open-source, and cross-platform
• Provides tools for automation and reproducibility
• Seamlessly integrates spatial data analysis and statistical analysis in one environment
• Handles vector and raster data

R Drawbacks for GIS Work

• R is less interactive than desktop applications like ArcGIS or QGIS
• Handling coordinate reference systems is more challenging
• In-memory analysis can be prohibitive for large data
• Steep learning curve

A Motivating Example

• Simple, expressive code

library(tmap)
library(tmaptools)
library(dplyr)

# find my town!
my_town <- tmaptools::geocode_OSM("Corvallis OR USA", as.sf = TRUE)
glimpse(my_town)

Rows: 1
Columns: 9
$ query   <chr> "Corvallis OR USA"
$ lat     <dbl> 44.56457
$ lon     <dbl> -123.262
$ lat_min <dbl> 44.51993
$ lat_max <dbl> 44.60724
$ lon_min <dbl> -123.3362
$ lon_max <dbl> -123.231
$ bbox    <POLYGON [°]> POLYGON ((-123.3362 44.5199...
$ point   <POINT [°]> POINT (-123.262 44.56457)

Your Turn

Our Turn

A Choropleth Map

• A choropleth map is a map that connects statistical summaries to geographic characteristics
• The tigris package can retrieve census, county, and state boundaries as vector sf data

library(tigris)
library(sf)
counties <- counties("Oregon", cb = TRUE)
counties$area <- as.numeric(st_area(counties))
glimpse(counties)

tm_shape(counties) +
  tm_polygons("area", style="quantile", title="Area of Oregon Counties")

Your Turn

library(mapview)
mapviewOptions(fgb=FALSE)
mapview(my_town, col="red", col.regions = "red") +
  mapview(counties, alpha.regions = .1)

Our Turn

Figure 2: Corvallis among distribution of Oregon counties.

Spatial Data Structures

• Spatial data structures are primarily organized through:

1. GDAL link here: For raster and feature abstraction and processing
2. PROJ link here: For coordinate transformations and projections
3. GEOS link here: For spatial operations (e.g., calculating buffers or centroids) on data with a projected coordinate reference system (CRS)

• We explore these core libraries throughout the workshop

Types of Spatial Data

• Vector data are comprised of points, lines, and polygons that represent discrete spatial entities (e.g., river, wtaershed, stream gauge)
• Raster data divide space into small rectangles (pixels) that represent spatially continuous phenomena like elevation or precipitation

Types of Spatial Data

Figure 3: Vector and raster data.

Vector Data Model

• Vector data are described using simple features
• Simple features is a standard way to specify how two-dimensional geometries are represented, stored in and retrieved from databases, and which geometric operations can be performed
• Provides framework for extending spatial POINTS to LINES and POLYGONS

Types of Spatial Data

Figure 4: The simple features data model.

The sf Package

• We use the sf package to work with spatial vector data in R
• Learn more about sf at their website r-spatial.github.io/sf and/or by running

vignette(package = "sf") # see which vignettes are available
vignette("sf1")          # open the introduction vignette

sf Primitives

• There are three basic geometric primitives with corresponding sf functions to create them:

1. POINT: created with st_points()
2. LINESTRING: created with st_linestring()
3. POLYGON: created with st_polygon()

• More complex simple feature geometries are built up using st_multipoint(), st_multilinestring(), st_multipolygon(), and st_geometrycollection()
• Building blocks for more complicated/large data structures

POINTs

• Points in sf are composed of coordinate pairs and a coordinate reference system
• Points have no length or area

# create sf object as point
pt1 <- st_point(c(1,3))
pt2 <- st_point(c(2,1))
pt3 <- st_point(c(3,2))
sf_point <- st_sfc(pt1, pt2, pt3)

ggplot(sf_point) + 
  geom_sf(color = "red") +
  labs(x = "X", y = "Y")

LINESTRINGs

• An ordered sequence of two or more POINTs
• Points in a line are called verticies (or nodes)
• Linestrings have length but no area

# create sf object as linestring
line1 <- sf::st_linestring(rbind(pt1, pt2))
line2 <- sf::st_linestring(rbind(pt2, pt3))
sf_linestring <- st_sfc(line1, line2)

ggplot(sf_linestring) + 
  geom_sf() +
  labs(x = "X", y = "Y")

POLYGONs

• Four or more points with the same starting and ending point
• Has length and area

pt4 <- pt1
poly1 <- st_polygon(list(rbind(pt1, pt2, pt3, pt4)))
sf_polygon <- st_sfc(poly1)

# dimension
sf::st_dimension(sf_polygon)

[1] 2

ggplot(sf_polygon) + 
  geom_sf(fill = "grey")

Putting Them All Together

ggplot() +
  geom_sf(data = sf_point, color = "red") +
  geom_sf(data = sf_linestring) +
  geom_sf(data = sf_polygon, fill = "grey")

Spatial Data Frames

• sf objects combine data.frame attributes with POINT, LINESTRING, and POLYGON building blocks (these are the geometry)
• Also characterize dimensionality, bounding box, coordinate reference system, attribute headers
• Read from shapefiles using st_read() or read_sf()

library(spData)
data("us_states")
print(us_states)

Spatial Data Frames

Figure 9: The simple features data structure.

Your Turn

Our Turn

plot(us_states["AREA"])

Raster Data Model

• Raster data can be continuous or categorical
• Can be image based; have a temporal component

Raster Data Model

Figure 11: The raster data model.

The terra Package

• We use the terra package to work with spatial raster data in R
• terra is the modern replacement to raster
• Learn more about terra at their website https://rspatial.org/
• stars for spatio-temporal rasters

Creating a Raster Object

• Create an empty SpatRaster object; define rows, columns, bounding box

library(terra)
myrast <- rast(ncol = 10, nrow = 10, xmax = -116, xmin = -126, ymin = 42, ymax = 46)
str(myrast)

S4 class 'SpatRaster' [package "terra"]

print(myrast)

class       : SpatRaster 
dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
resolution  : 1, 0.4  (x, y)
extent      : -126, -116, 42, 46  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 

Your Turn

Our Turn

class(myrast)

[1] "SpatRaster"
attr(,"package")
[1] "terra"

ncell(myrast)

[1] 100

res(myrast)

[1] 1.0 0.4

ext(myrast)

SpatExtent : -126, -116, 42, 46 (xmin, xmax, ymin, ymax)

crs(myrast)

[1] "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]],\n        ID[\"EPSG\",6326]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433],\n        ID[\"EPSG\",8901]],\n    CS[ellipsoidal,2],\n        AXIS[\"longitude\",east,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433,\n                ID[\"EPSG\",9122]]],\n        AXIS[\"latitude\",north,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433,\n                ID[\"EPSG\",9122]]]]"

ncell(myrast)

[1] 100

Manipulating Raster Objects

• Let’s add values to the raster object

values(myrast) <- 1:ncell(myrast)

plot(myrast)

Reading a Raster Object

• Lets read in a raster object from the spDataLarge package that contains elevation data from Zion National Park

raster_filepath <- system.file("raster/srtm.tif", package = "spDataLarge")
my_rast <- rast(raster_filepath)

plot(my_rast)

Your Turn

landsat <- system.file("raster/landsat.tif", package = "spDataLarge")
landsat <- rast(landsat)

Our Turn

plot(landsat)

Figure 14: Raster plot of Landsat 8 data from Zion National Park.

Spatio-Temporal Raster Objects

• The stars package has tools for spatio-temporal raster data
• stars integrates with sf (recall that terra does not)
• Learn more at their website https://r-spatial.github.io/stars/

Coordinate Reference Systems

A coordinate reference system (CRS) is made up of several components

1. Coordinate system: The x,y grid that defines where your data lies in space
2. Horizontal and vertical units: The units that describe the coordinate system
3. Datum: The modeled version of the earth’s shape
4. Projection details (if applicable): The mathematical equation used to flatten objects from round surface (earth) to flat surface (paper or screen)

Properly specifying coordinate reference systems is a crucial step in handling or analyzing spatial data

The Ellipsoid and Geoid

• The earth can be thought of as an ellipsoid defined by a semi-major (equatorial) axis and a semi-minor (polar) axis
• More precisely, it is a geoid that is not perfectly smooth
• A datum is built on top of the ellipsoid and incorporates local features
• The most appropriate datum depends on the location of interest (in the US we typically use NAD83)

The Ellipsoid and Geoid

Figure 15: The ellipsoid and geoid.

Your Turn

Specifying a CRS

There are many ways to specify a CRS:

1. Proj4
   • +proj = longlat + ellps=WGS84 ...
2. OGC WKT / ESRI WKT
   • GEOGCS["WGS 84", DATUM["WGS_1984, SPHERIOD["WGS 84 ...]]]
3. authority:code identifier (modern/preferred)
   • EPSG: 4326

Projected Coordinate Systems

• Projected coordinates have been projected to two-dimensional space according to a CRS
• Projected coordinates have an origin, x-axis, y-axis, and unit of measure
• Conformal projections preserve shape
• Equal-area projections preserve area
• Equidistant projections preserve distance
• Azimuthal projection preserves direction

Projected Coordinate Systems

• Here an example using vector data; see the workbook for an example using raster data

library(Rspatialworkshop)
data("pnw")
# transform one to the other
utmz11 <- 2153
pnw <- st_transform(pnw, crs = utmz11)
ggplot() + 
  geom_sf(data = pnw,  color="black", fill=NA) +
  labs(title="State Boundaries in the Pacific Northwest") +
  theme_bw() 

Your Turn

Our Turn

pnw2 <- st_transform(pnw, crs = 5070)
plot(st_geometry(pnw2))

A Note on S2

• sf version 1.0.0 supports spherical geometry operations via its interface to Google’s S2 spherical geometry engine
• S2 is an example of a Discrete Global Grid System (DGGS)
• sf can run with s2 on or off and by default the S2 geometry engine is turned on:

sf::sf_use_s2()

Geographic Data I/O (Input/Ouptut)

There are several ways to read spatial data into R

1. Load spatial data from our machine or a remote source
2. Load spatial data as part of an R package
3. Load data using an API (which often makes use of an R package)
4. Convert flat files with x, y coordinates to spatial data
5. Geocoding data “by-hand” (we saw this earlier)

Vector Data I/O

sf can read numerous file types:

1. shapefiles
2. geodatabases
3. geopackages
4. geojson
5. spatial database files

Read in a Shapefile

filepath <- system.file("extdata/city_limits.shp", package = "Rspatialworkshop")
citylims <- read_sf(filepath)
plot(st_geometry(citylims), axes = T, main = "Oregon City Limits")

Your Turn

Our Turn

Read in a Geodatabase

filepath <- system.file("extdata/StateParkBoundaries.gdb", package = "Rspatialworkshop")
# List all feature classes in a file geodatabase
st_layers(filepath)
# Read the feature class
parks <- st_read(dsn = filepath, layer="StateParkBoundaries")
ggplot(parks) + 
  geom_sf()

Read in a Geopackage

filepath <- system.file("gpkg/nc.gpkg", package="sf")
nc <- read_sf(filepath)
ggplot(nc) + 
  geom_sf()

Geopackage Advantages

Why use geopackages over shapefiles?

1. Geopackages avoid mult-file format of shapefiles
2. Geopackages avoid the 2gb limit of shapefiles
3. Geopackages are open-source and follow OGC standards
4. Geopackages are lighter in file size than shapefiles
5. Geopackages avoid the 10-character limit to column headers in shapefile attribute tables (stored in archaic .dbf files)

Other Approaches

The workbook shows how to read in data using

1. Open spatial data sources
2. R packages
3. OpenStreetMap

Raster Data I/O

• Read in data using rast()

filepath <- system.file("ex/elev.tif", package="terra")
elev <- rast(filepath)
barplot(elev, digits = -1, las = 2, ylab = "Frequency")

Raster Data I/O

plot(elev)

Figure 21: Elevation Data

Raster Data I/O

library(stars)
filepath <- system.file("tif/L7_ETMs.tif", package = "stars")
L7 <- read_stars(filepath) |>
  slice(index = 1, along = "band")
plot(L7)

Convert Flat Files to Spatial

It is common for a flat file (e.g., a .csv) to have columns that indicate coordinates – how do we turn this into a spatial vector object?

filepath <- system.file("extdata/Gages_flowdata.csv", package = "Rspatialworkshop")
gages <- read.csv(filepath)
gages_sf <- gages |>
  st_as_sf(coords = c("LON_SITE", "LAT_SITE"), crs = 4269, remove = FALSE)
ggplot(data = gages_sf) + 
  geom_sf()